// STL
#include <cmath>                // std::sqrt()
#include <stdexcept>            // std::runtime_error()
#include <utility>              // std::pair<>, ::make_pair()

// Boost
#include <boost/math/distributions/normal.hpp> // boost::math::normal

// stats++
#include "statsxx/sampling.hpp" // statsxx::sampling::sample_size_m()


//
// DESC: Determine the significance level for simultaneously estimating the parameters of a multinomial population within a distance of the true values at a sample size.
//
// -----
//
// INPUT:
//
//     d^2n  : d == distance, n == sample size
//
// -----
//
// OUTPUT:
//
//     alpha : significance level
//     m     : the number of nonzero parameters, at which the "worst case" occurs for the given alpha
//
// -----
//
// NOTE: See:
//
//     S. K. Thompson, "Sample Size for Estimating Multinomial Proportions," The American Statistician 41, 42--46 (1987)
//
// -----
//
inline std::pair<
                 double, // alpha
                 int     // m
                 > statsxx::sampling::significance_multinomial(
                                                               const double d2n
                                                               )
{
    double alpha;
    int    m;

    // INITIALIZATION
    boost::math::normal dist;

    // SELF CONSISTENTLY FIND m
    m = 6;

    do
    {
        // FOR THIS m, CALCULATE z
        double z = d2n*m/(1. - 1./m);
        z = std::sqrt(z);

        // CALCULATE alpha FROM z
        alpha = (2*m)*(1. - boost::math::cdf(dist, z));

        // NOTE: In the case of insufficient data and/or small d, alpha can sometimes be larger than 1.
        if( alpha >= 1. )
        {
            // NOTE: ... So we return (correctly), 1 alpha for zero categories.
            return std::make_pair(
                                  1.,
                                  0
                                  );
        }
        // NOTE: ... And in the case of significant data and/or large d, alpha can sometimes be 0.
        else if( alpha == 0. )
        {
            // NOTE: ... So we return (correctly), 0 alpha for at least 1 category.
            return std::make_pair(
                                  0.,
                                  1
                                  );
        }



        // CHECK FOR SELF CONSISTENCY
        if( statsxx::sampling::sample_size_m(alpha) != m )
        {
            m--;

            if( m == 1 )
            {
                throw std::runtime_error("error in statsxx::sampling::significance_multinomial(): no self-consistent solution found");
            }
        }
        else
        {
            break;
        }

    } while( true );

    // RETURN
    return std::make_pair(
                          alpha,
                          m
                          );
}
